# konrad.smolinski@gmail.com
# date:		   23/11/2010
# last update: 1/12/2010
#
# info:
#	we assume:
#	 Y is binary: Y in {0,1}
#	 X has K-points of support: {x_1,...,x_K}
# 	 Z has R-points of support: {z_1,...,z_R}
# 	 U = (U1,U2)
#
# dependencies:
# library: gss
# sources: rcWatchOut.R
# dependencies: rcNormalCopula()
# ---------------------------------------------------------------------------------------
# rcBaseIntegrate() : integrate Normal Copula over the Base Components 
rcBaseIntegrate <- function(a=c(-0.1,0.5),b=c(0.1,1),xval=c(-0.5,0.5),mu=c(0,0),quadPt=smolyak.quad(2,35)){
	# parameterization
	sigMat <- cbind(c(1,b[1]),c(b[1],b[2]+b[1]^2))	
	
	# quadratures:	
	pts <- quadPt$pt; ww <- quadPt$wt
	
	# Integrand evaluated at sparse grid points
	val <- NULL
	for(i in 1:nrow(pts)) val <- c(val,rcNormalCopula(pts[i,],mu=mu,sig=sigMat) ) 
	int <- val*ww  

	# important Part:
	# support and pie slices in R^2
	ysupp <- c(FALSE,TRUE)
	K <- length(xval)

	# Y=0 conditional on X=x[k]
	xy0 <- NULL
	for(i in 1:K) xy0 <- cbind(xy0, (pts[,2] < pnorm(a[1]+(a[2]+qnorm(pts[,1],sd=sqrt(b[2]+b[1]^2)))*xval[i],mean=mu[2],sd=1) ))

	# ----
	baseCompPtsMat <- NULL
	xysupp <- list()
	for(i in 1:K) xysupp[[i]] <- ysupp

	baseComponents <- rcBase(K)
	nOfBaseComponents <- nrow(baseComponents)

	# this can be done much shorter and more neat 
	for(i in 1:nOfBaseComponents){

		if(baseComponents[i,1]==0) baseCompPts <- xy0[,1] 
		if(baseComponents[i,1]==1) baseCompPts <- !xy0[,1] 

		for(j in 2:K){
			if(baseComponents[i,j]==0) temp <- xy0[,j] 
			if(baseComponents[i,j]==1) temp <- !xy0[,j] 
			baseCompPts <- (baseCompPts & temp)
		}
		baseCompPtsMat <- cbind(baseCompPtsMat,baseCompPts)
	}
	
	# Bind integrated Values:
	intBaseComponents <- NULL
	for(i in 1:nOfBaseComponents)	intBaseComponents <- rbind(intBaseComponents,sum(int[baseCompPtsMat[,i]]))

return(intBaseComponents)
}
# --------------------------------------------------------------------------
# Validate:
codeVal <- FALSE
if(codeVal){
cat("\n Validate: rcBaseIntegrate() : \n\n")
intVec <- rcBaseIntegrate()
intVec
# Should sum to 1:
sum(intVec)
}
